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When considering NLO corrections to thermal particle production in the "relativistic" regime, 
in which the invariant mass squared of the produced particle is /C 2 ~ (ttT) 2 , then the produc- 



tion rate can be expressed as a sum of a few universal "master" spectral functions. Taking 
the most complicated 2-loop master as an example, a general strategy for obtaining a conver- 

CN , gent 2-dimensional integral representation is suggested. The analysis applies both to bosonic 

and fermionic statistics, and shows that for this master the non-relativistic approximation is 

^^ , only accurate for /C 2 >(8irT) 2 , whereas the zero-momentum approximation works surprisingly 

well. Once the simpler masters have been similarly resolved, NLO results for quantities such 
as the right-handed neutrino production rate from a Standard Model plasma or the dilepton 
production rate from a QCD plasma can be assembled for /C 2 ~ (ttT) 2 . 
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1. Introduction 

Given the remarkable confirmation of microscopic Standard Model and QCD physics through 
the Large Hadron Collider program, it appears well-motivated to work out how the same 
interactions behaved in the macroscopic thermal environment of the Early Universe. Whereas 
there are good tools available for doing this for thermodynamic quantities such as the overall 
equation of state, much less is known about real-time rates, such as particle production rates 
or the equilibration rates related to the most weakly interacting degrees of freedom. Indeed, in 
most cases only phenomenological estimates or leading-order (LO) weak-coupling expressions 
are available. Given the known infrared problems of thermal field theory, which imply that 
next-to-leading order (NLO) corrections may be surprisingly large, further work is needed in 
order to get an impression on the numerical accuracy of the results currently available. 

In relativistic thermal field theory, the structure of the weak-coupling expansion depends 
sensitively on the physical scales of the problem. Focussing on the simplest situation, in 
which there is only one zero-temperature mass scale, denoted by the invariant mass squared 
K? > 0, there are different "regimes" depending on the ratio K? /(ttT) 2 . 1 In the so-called 
non-relativistic regime, K? 3> (irT) 2 , results can be represented in the form of an Operator 
Product Expansion [1], and thermal corrections are in general small (power-suppressed). In 
contrast, in the so-called ultrarelativistic regime, K? <C (itT) 2 , the naive loop expansion 
breaks down, and extensive resummations (incorporating for instance the physics of the 
Landau- Pomeranchuk-Migdal effect) are needed for determining even the LO result. The 
current technology for this has been developed in the context of the photon production rate 
from a QCD plasma [2, 3, 4] 2 and subsequently reformulated and applied to the right-handed 
neutrino production rate from a Standard Model plasma [8, 9, 10]. (Very recently the photon 
production rate has been determined up to NLO [11], which here means 0(a 1 ' 2 ) rather than 
the usual O(a), where a = <7 2 /(47r).) 

The focus of the present paper is the so-called relativistic regime, K? ~ i^T) 2 . Such 
results may permit for an interpolation between the non-relativistic and ultrarelativistic cases, 
thereby perhaps yielding phenomenologically broadly applicable expressions (the cosmological 
evolution is illustrated in fig. 1). At NLO, which in the relativistic regime means 0(a), no 
infrared divergences are expected to be encountered in the full result, and there is no need 
for resummations either. However, because of the loss of Lorentz symmetry within a thermal 
bath, the analysis is technically cumbersome. This is the case particularly if a non-zero 
momentum k = |k| ^ is considered, in which case even LO expressions are non-trivial [12]. 
Nevertheless, as will be demonstrated, NLO results can be worked out with some effort. 

Examples of concrete applications that we have in mind are the right-handed neutrino 



x To be precise, because of the absence of Lorentz invariance within a heat bath the magnitudes of k± 
(ko ± k)/2 need to be separately compared with nT, and the number of regimes may proliferate. 
2 Similar resummations are needed for dilepton production in the soft regime K? <C (tvT) 2 [5, 6, 7]. 
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Figure 1: An example of how during the cosmological evolution, different kinematic regimes are 
crossed in the thermal production of a particle species with a fixed invariant mass K? = M 2 . A 
generic weak gauge coupling is denoted by g [a = g 2 /(in)], and fc av is defined by eq. (8.3). If a fixed 
invariant energy (E + E ,) 2 — (k + k ,) 2 ~ I GeV 2 of a dilepton pair is considered, a similar 
evolution takes place during a heavy ion collision, however the temperature range is rather narrow. 



production rate in the Early Universe, motivated by leptogenesis or dark matter computa- 
tions [13, 14], and the dilepton production rate in hot QCD, motivated by heavy ion collision 
experiments. For the former, the NLO level has already been reached in the non-relativistic 
regime [15, 16]. For the latter, the NLO level was reached long ago at vanishing momentum 
in the relativistic regime [17, 18, 19]. More recently, the case of a vanishing momentum in 
the relativistic regime has been studied more generally [20], expressing results in terms of 
a finite number of universal "master" structures; in particle language the cases considered 
were the (off-shell) thermal axion and dilaton [20] or graviton [21] production rates. (The 
physical relevance of these computations is related to lattice determinations of the corre- 
sponding "transport coefficients" , which are the rate of anomalous chirality violation and the 
bulk and shear viscosities. Indeed, because of the necessity of analytic continuation from Eu- 
clidean signature, ultraviolet features need to subtracted from non-perturbative data before 
an extrapolation to ko — > becomes possible even in principle [22, 23, 24].) 

The purpose of the present paper is to evaluate the spectral function corresponding to the 
most complicated 2-loop master topology (fig. 2) at non-zero momentum in the relativistic 
regime. This topology has merited extensive investigations under other circumstances. In- 
deed, at zero temperature in 3 space-time dimensions even the massive case can be solved [25], 
whereas in the physical 4-dimensional case the result has an extremely rich structure [26] 



which, despite a vast body of work, still remains under further investigation today (see, e.g., 
ref. [27] and references therein). In the massless limit, however, the result vanishes at zero 
temperature. Nevertheless the spectral function (cut) is non-zero and possesses a rich struc- 
ture at finite temperature. Previous thermal analyses exist at zero momentum [20] as well as 
at zero energy but non-zero momentum [28], the latter case corresponding to the space-like 
domain relevant for the physics of plasma screening. 

Recently, a separate line of study of thermal particle production in the relativistic regime 
has been initiated whose goals appear to partly overlap with those of the present paper [29]. 
The authors give a 4-dimensional integral representation for the gauge boson contribution to 
the right-handed neutrino production rate, but no numerical evaluation is shown. Given the 
technical complexity of the problem, it appears welcome that two independent and method- 
ologically different computations are being pursued, permitting in the end to crosscheck the 
validity of both results. 

The plan of this paper is the following. We start by recalling how NLO results for two 
physical observables can be expressed in terms of a finite number of simple master sum- 
integrals (sec. 2). After carrying out Matsubara sums, the cut of the most complicated 
master can furthermore be decomposed into processes representing real and virtual corrections 
(sec. 3). The real corrections are analyzed in sec. 4, showing that through a suitable choice 
of variables they can be given a 2-dimensional integral representation (soft and collinear 
divergences are regulated by an auxiliary mass parameter at this stage); the same task is 
accomplished for the virtual corrections in sec. 5. Both real and virtual corrections are 
divergent if the auxiliary mass parameter is sent to zero; in sec. 6 it is shown that the sum 
remains finite. Finally a form suitable for practical evaluation is given in sec. 7, whereas 
numerical comparisons with known limiting values comprise sec. 8. A brief summary and 
outlook are offered in sec. 9. 



2. Physical observables 

To underline the significance of the master spectral function considered, we start by recalling 
two separate physics contexts in which it plays a role. 

To leading order in a em = e 2 /(4ir), the production rate of fx~ fi + pairs from a hot QCD 
medium can be expressed as [30, 31, 32] 
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where Q q is the electric charge of quark of flavour q in units of e, n B is the Bose distribution, 
and 



Pqq ,{K) = j^ e ilcx (± [JfWJq'M ) (2-2) 



T 
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is the spectral function corresponding to the vector current. (The metric convention /C ■ X = 
kQXQ — k • x is assumed.) The spectral function can in turn be expressed as the cut of the 
corresponding Euclidean correlator, 



Pqq'(fc) = Im U Eqq/ (K) 



k n — >— i[ko+iQ + ] 



(2.3) 



Here K = (fc ra ,k) denotes a Matsubara four-momentum, with K 2 = k 2 + k 2 , k = |k|. 

The Euclidean correlator can be computed with standard path integral techniques. After 
carrying out the Dirac algebra and making use of substitutions of integration variables, it 
can be expressed as a sum of a finite number of "master" structures. Up to 2-loop level, the 
expression reads 
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Here D is the dimensionality of space-time, and C P = (N 2 — i)/(2N c ). According to standard 
conventions, J{ P} denotes a sum-integral with fermionic Matsubara momenta, and % P one 
with bosonic Matsubara momenta. (The generalization of eq. (2.4) to a finite quark mass 
can be found in ref. [33].) 

The different terms of eq. (2.4) can be referred to as master sum-integrals; their cuts, defined 
in accordance with eq. (2.3), are the corresponding master spectral functions. In particular, 
the most complicated master sum-integral (the only one with 5 different propagators) is 
defined from now on as 

H K ) = j™^ ^™,n m5 I K lr, w, n m , (2-5) 
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and is illustrated graphically in fig. 2. Note that a mass regulator A has been introduced for 
intermediate stages of the computation. The subscript (...)j corresponds to the labelling of 
the various masters in ref. [16]. The statistics of P,Q are meant to be interpreted openly, 
and are labelled by the indices aQ...a^ as shown in fig. 2. 
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Figure 2: The master Feynman diagram considered, denned by eq. (2.5). The integration measure 
corresponds to a Matsubara sum-integral, and the spectral function is the discontinuity (imaginary 
part) of the result once the external Matsubara frequency is analytically continued as k n — > — i[ko + 
i0 + ]. Dashed lines represent massless propagators; the wiggly line is regulated through a small mass 
parameter, A, for intermediate stages of the computation. The internal and external lines could be 
either fermions or bosons [denoted by a — —1 and +1, respectively, with fi £ {0, ...,5}], however 
fermion number conservation is assumed at each vertex. 

The indices aQ...a^ take the value +1 for bosons and —1 for fermions. Assuming fermion 
number conservation at each vertex, only three of the indices are independent; we can choose 
<7i,o"4,0"5 to play this role. Then 

(To = 0"1<74 i Ol = 0"lO"5 i <73 = C4C5 ■ (2-6) 

The case appearing in eq. (2.4) corresponds to (a"iU40" 5 ) =( h), and a numerical evaluation 

for the corresponding spectral function, defined in eq. (3.1), is shown in fig. 8 below [along 
the curve in the (k, A:o)-plane shown in fig. 1]. 

A completely different physics application, leading to the same master spectral functions 
but with (partly) different statistics, is that of right-handed neutrino production in the Early 
Universe. Again, a gauge-invariant Euclidean correlator can be defined, 

n E (*r) = Tr [i$ [ I'' T dr f e* K - x ((ftaJ)(X) (iaj)(0))] } , (2.7) 

where X = (t, x); I is a lepton doublet; <fi is a Higgs doublet; and a h , a R are chiral projectors. 
The production rate now reads 



dAU/C) 2 \h, J 2 



'R 



d 4 ^d 3 k (2TT) 3 k 



n F (k )p(fC) , (2i 



where n F is the Fermi distribution; h u a bare neutrino Yukawa coupling; and ko = \Jk 2 + M 2 , 
with M denoting the right-handed neutrino mass. The spectral function p is obtained from 
the cut of the Euclidean correlator IT E just like in eq. (2.3). 



The 2-loop expression for II E can be written in a form analogous to eq. (2.4) [16]: 
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(2.9) 



Here A^ is the Higgs self-coupling, ht is the top Yukawa coupling, and g 1 , g 2 are the hyper- 
charge and weak gauge coupling, respectively. 

The master sum- integrals appearing in eq. (2.9) are identical to those in eq. (2.4), apart 
from their statistics. Noting that K is now fermionic, the most complicated case, eq. (2.5), 

appears with the indices (aiders) = (-1 h). A numerical evaluation is presented in fig. 8 

along the curve in the (/c, feo)-plane shown in fig. 1. 

The two examples discussed should serve as illustrations of concrete applications of the 
spectral function related to eq. (2.5), but do not exhaust the cases considered in the literature. 
In particular, physical observables leading to the case (aiders) = (+ + +) have been analyzed 
in some detail for k = in the relativistic regime in refs. [20, 21], and for k ^ in the non- 
relativistic regime in refs. [34, 35, 36]. 

As a final remark we mention that above the physical observable was a Lorentz scalar, and 
subsequently the master structures are scalars as well. In the right-handed neutrino case it 
may be of interest to ultimately compute the whole self-energy matrix. In this case tensor 



sum-integrals appear apart from scalar ones; the corresponding results have been worked out 
in the non-relativistic regime in ref. [37] . A general discussion concerning the thermal tensor 
basis can be found in ref. [38]. 



3. Splitup into real and virtual corrections 

Given eq. (2.5), the corresponding spectral function can be defined as 

Pt. = Im { X j}fc, l ^-*[fc ( ,+*o+] • (3.1) 

It is well established at 1-loop level that the cut yields a structure reminiscent of a Boltzmann 
equation, with scattering amplitudes squared multiplied by appropriate phase space distri- 
butions [39]. A similar result applies at 2-loop level, but is somewhat complicated by the 
fact that now virtual particles also appear. In any case, making use of standard techniques, 
explained in some detail for instance in appendix A of ref. [40], and shifting four-momenta, 
the following result can be obtained: 
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Here P refers to principal value integration which renders changes of integration variables 
unproblematic; and 
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Furthermore the phase space distributions are defined as 
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The notation (...) |o =±e implies that terms with both signs be summed together. The chan- 
nels labelled by (rl)-(r4) are referred to as real corrections] (vl)-(v3) as virtual corrections. 
The corresponding physical processes are illustrated in fig. 3. 

The last line of eq. (3.2) adds terms with ko — > —ko, which can be shown to render the 
spectral function antisymmetric in this exchange. In the following we restrict to ko > 0, 
and in this case only the terms (rl)— (r4) and (vl)-(v3) contribute. (The unlabelled virtual 
corrections below (v3) in eq. (3.2) only contribute in the space-like domain.) 



4. Real corrections 

Many numerical evaluations of phase space integrals like (rl)-(r4) can be found in the liter- 
ature, but we are not aware of a previous representation as a 2-dimensional integral. As is 
now demonstrated, a fairly explicit expression can be obtained even for A 7^ 0. 
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Figure 3: Graphical illustrations of the cut of fig. 2, given explicitly in eq. (3.2). Only the decay 
channels are shown, but as the phase space distributions in eq. (3.2) indicate, the inverse processes are 
included as well. The processes (r2), (r3), (r4) represent interference terms between (u,s), (t,s) and 
(t, u) channels. For the virtual corrections the labelling (vl)-(v3) refers to cases in which different 
internal lines of the closed loop have a thermal weight. 



4.1. Integration variables 



As is common in particle kinematics, a fruitful approach is to represent a complicated phase 
space as a convolution of simpler ones [41]. Taking case (rl) from eq. (3.2) as an example, 
we may rewrite it as 
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where in the first step we substituted p — > k — p and integrated over r; and in the second 
step introduced a variable po implementing the convolution. The function $ rl contains phase 
space distributions as well as propagators; the symbol "•" stands for variables not shown 
explicitly. 

The key observation now is that the first three arguments of <l> rl , which appear in the phase 
space distributions, do not contain the variable p. Therefore, it is advantageous to take po 
and q as the outer integration variables; and p as well as one azimuthal angle that is not 
fixed by the 5-constraints as the inner ones. The inner integrations contain no phase space 
distributions and, as we will see, can be carried out explicitly. They also turn out to capture 
collinear phase space singularities in a manageable form. 

The only challenge with this strategy is that it is tedious to work out the limits of the p- 
integration in the plane (po, q). Nevertheless, with some work, this challenge is surmountable. 
Employing the labelling shown in fig. 4, the ranges for the case (rl) are [k± = (ko ± k)/2]: 
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Figure 4: Integration ranges in the (po,<z)-plane f° r rea l corrections (for A = fco/20). The delimiting 
curves are q = k + - A 2 /[4(p - fe+)], q = k- - A 2 /[4(p - fe-)], # = Po — k- - A 2 /(4fc_), g = 
Po — k + — A 2 /(4fc + ), as well as q = po ± A. The various channels, separated by the dotted lines and 
the gap, are in the middle (rl); at po > k (r2); at q < (r3); and at q > p (r4). 



As it turns out, the situations ko < 2k + \/k 2 + A 2 and ko > 2k + Vk 2 + A 2 need to be handled 
separately, and this leads to the two cases (c) and (c). 
With the ranges at hand, the integrals over the angles 



cosx 



pk 



cos# 



pq 



pk pq 

can be carried out in eq. (4.1), thereby removing the two (^-functions. We are left with 



(4.10) 



1 



/■Pmax(po,<?) 

/Vi "= 77 N 3fc / d Po dq dp ( $ Tl (k - p \q\po ~ q\ 



(4.11) 



Here (...) refers to an azimuthal average, and O rl is composed of the domains (a)-(e) of fig. 4. 
It is a nice crosscheck to set <3? rl —> 1 and carry out the remaining integrals. In the 
absence of phase space distributions, the integral is Lorentz-invariant and can alternatively 
be performed in a suitable frame to find a simple result: 



Pril$ rl _».i 



(2^) A 5^5{-K + V + Q + 1Z) 1 



p,q.r 



8e p e q E r 



(4tt) 



K, 



A 4 A 2 , A 2 

1 In — 

4 4/C 2 2 K 2 



(4.12) 



This can be reproduced from eq. (4.11), both for ko < 2k + \Jk 2 + A 2 and ko > 2k + \Jk 2 + A 2 . 
However, in the thermal case Lorentz symmetry is not particularly helpful because the plasma 
defines a special frame, and we need to make use of eqs. (4.2)— (4.9). 
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The other channels can be handled similarly. For the case (r2), the shift p — > p — k and a 
subsequent integration over r yields 



(2ir) 4 5(-k -e p + e q + E r )5^(-k - p + q+ r) 

8e„eK ~r2V«=pl*,l"rl 



Pr2 = / — £ £— £ — % 2 {e v \e q \E 

'p,q,r ot p t q I1/ r 

(■OO 



f 2ir5(-ko - e k + po)5{po - e q - E ) . , . . ., 

dpo / t ® l2 {p - k \e q \p - e q \-) 

oo Jp,q ° e pfc e g- c 'gp 






/• rPmax(P0,<?) / . 

/ dp dq dp ($ l2 (po - ko\q\po - q\-)\ . (4.13) 

./fi,o Jpmln(PO,Q) 



The domain is displayed in fig. 4, with the ranges 



(f): pe (2fc + -po,Po-2fc_J , (4.14) 

(g): p€ (g-V(9-Po) 2 -A 2 ,p -2fc_) , (4.15) 



(h): p€ po-2fc+,po-2fc- , (4.16) 



(i): pe (- (7 + V(9-Po) 2 -A 2 ,p -2A ; _j . (4.17) 

For the case (r3), the shifts p — > k — p, q — > — q and a subsequent integration over r yield 
f (27r) 4 <5(-A;o + 6 p -e g + J B r ) ( 5( 3 )(-k + p-q + r) 

Pr3 = / " J— p ^(^M^rlO 

» /• 27T(5(-fco + e fc +po)(J(po + e,-£? ffl ,) 

dpo / ^ = *r3("o-Po|e,|Po + e,|-) 

oo ^p,q ° e pk e q qp 



\ f rPmax(P0,9) 

(4tt) 3 A; .,<>,, ./ /)|!:;!lU , iM/i 



r rPmax(po,g) , . 

/ dp dq dp ($ l3 (k - Po\ -q\po-q\-)) ■ (4.18) 

Jn r3 Jpmln(PO,q) 



In the 2nd step we substituted formally e„ — > — q, which permits us to represent the domain 
as displayed in fig. 4, with the ranges 

GQ: P e (« + V(9 - Po) 2 - A 2 , -g + y/(q - p ) 2 - A 2 ) , (4.19) 



(g) : p£ [2k- - po, -q + V(q-Po) 2 -^) , (4-20) 

(h): p£ (2k.-p ,2k + - Po ) , (4.21) 



(i) : p G (g + V(9-Po) 2 -A 2 , 2A;+ - Po j • (4.22) 

12 



Finally, for the case (r4) , the shift p — > k — p and a subsequent integration over r yields 

Pr4 = / S77~? $T4( e p\ e q \ E r\-) 

./p,q,r otp6g-C/ r 



, 2vrJ(-/co + e k + Po )S(po - e + £7 ) 

d Po / ^ ^ ^4(^0 -Po|e |e - Pol 



'" ' - — — p *r4V"'U -^U|c g |c g 

00 Jp,q "^pk^q^qp 



I I- / , Pmax(po,<?) / . 

= f 4 N 3fc / d Podg/ dp^ r4 (/c -Pokk-Pol-)/ • ( 4 - 23 ) 
The domain is displayed in fig. 4, with the ranges 

(j): pe (2/c_-p ,2/c+-po) , (4-24) 

(k): P£ ((?-V(9-Po) 2 -A 2 ,2A; + -po), (4.25) 



(k): pe 2fc_-po,<Z+V(9-Po) 2 -A 2 ) , (4.26) 



(1): p € g - y/(q - p ) 2 - A 2 , q + V(<7 - p ) 2 - A 2 . (4.27) 



4.2. Crossing symmetry 

At zero temperature, the four channels (rl)— (r4) of fig. 3 are related by a crossing symmetry, 
and it is interesting to see how the presence of phase space distributions in eq. (3.2) changes 
the situation. Inserting the arguments of <3? as displayed in eqs. (4.11), (4.13), (4.18), (4.23) 
into expressions obtained from eq. (3.2), we find 



$ Tl (ko - po\q\p Q - q\-)) - -{$ T2 (Po-ko\q\po-q\- 

$r3( fe o -Pol -q\po-q\- 
$ r4 (/co -Po\q\q~Po\ 
= n~^(ko) n a4 (k - p ) n ff2 (q) n a5 (p - q) 

xpf -j -V - \= *-\ 1, (4.28) 

I 2(Po ~ P 2 ) \ ( k ~ q) 2 ~ Zq k I (fe -P0|9|P0-?|-) J 

where on the last line the arguments (...|...|...|-) refer to e , , e q , and E qp , respectively, and 
we made use of the fact that the dependence on e , and E is quadratic [cf. eqs. (4.35), 
(4.36), (4.38) below] and that the sign of q inside e k plays no role. To arrive at eq. (4.28) the 
conservation of fermion number [cf. eq. (2.6)], leading to 02O5 = G\ etc, as well as an identity 
following from eq. (3.5), 

ve e ' T n ff (e) = -n ff {-e) , (4.29) 

were assumed. The universal form in eq. (4.28) implies that the azimuthal average (...) and 
the subsequent integration over p only need to be carried out for one single function. 
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4.3. Inner integrations 

We now consider the integrations still to be performed, i.e. 

/•Pmax(po,g) f AC 4 / 1 \ 1 

/(Po ' 9) "/ , , dpF 2fa 2 rMffcn ^ e 2 / • (430) 

• y P mi n(P0,<?) I 4P P ) \ l fc - q) t qk / (k -p \q\p -q\-) > 

Let us start with the azimuthal average. Parametrizing 3 

P = P (0,0,1) , (4.31) 

k = fc(sinx,0,cosx) , (4.32) 

q = q (sin 6 cos (p, sin 6 sin 99, cos 6*) , (4.33) 

the integral over 93 is readily carried out: 

f fepf * )=Ref: ign(a M , (4.34) 

where 

a = (fco — 5) —A; —g + 2/cq cos % cos , /3 = 2kq sin x sin . (4.35) 

Furthermore, according to eq. (4.10), the angles can be written as 

n 2 _i_ u2 2 

cosx= — ; — = 1H z— ; (4.36) 



2pA; 



p 2 + c? 2 + A 2 - £ 2 p 





(p- 


Po + 2fc_)(p + p - 


2k + ) 




1 

1 , (P 


2p/c 


-2k. 


-) 


+ 
1 


1 
(p- 


2p/c 
9) 2 -[(9-po) 2 -^ 


< 2 \ 




. (P 


2pg 
+ q) 2 - [(q - p ) 2 - 


■A 2 ] 





2pq 



(4.37) 



cosfl = - ^ = 1 + ^— ^ L ^ ^ uy lli (4.38) 

2gp 



(4.39) 



where e . , -E™, were inserted from the arguments shown in eq. (4.30). It can seen that at 
the boundaries of the p-integration, cf. eqs. (4.2)-(4.9), (4.14)-(4.17), (4.19)-(4.22), (4.24)- 
(4.27), one of the cosines evaluates to ±1, and correspondingly one of the sines vanishes. 
Therefore, at the boundaries the function j3 given in eq. (4.35) vanishes; this observation will 
turn out to be useful in a moment. 

Inspecting the expressions it is now possible to realize that, in general, the dependence of 
a 2 — /3 2 appearing in eq. (4.34) on p is of the form 

a 2 -(3 2 = ap 2 + b + ^ . (4.40) 

p2 



3 Note that this parametrization can be used both for positive and negative q. 
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Here a = (ko — q) 2 and the other coefficients are more complicated. Remarkably, this func- 
tional form implies that the integral defined in eq. (4.30) can be carried out: 



/C 4 fw 2 pdp f 

1{pthq) = Tj p . pJ^H 



sign (a J 



\/ap 4 + bp 2 + 



/C 4 sign (a) 



V°Po + & Po + 



In 



Pi - pL 



pI - pLk 



+ ln 



(V^O + bPo + C + \/ a P m a X + hPmax + rf ~ «(P0 ~ P 



2 ^2 

max / 



Wapt + bpl + c+ ^/apt~Tbpi~Tc) 2 - a(p 2 Q - p^J 2 
Furthermore, both square roots can be simplified: the prefactor contains the function 



(4.41) 



F(po,q) = yapo + bpl + c 

= V(9 " Po) 2 ^ 4 - A 2 [/C 2 + 4qpo - 2k (g + Po )]JC 2 + A 4 fc 2 , 



(4.42) 



whereas, as already mentioned in connection with eqs. (4.36)-(4.39), the function j3 vanishes 
for p = p min and p = p max . Therefore, according to eq. (4.40), 



V a Pmax + "Pmax + C — Pmax | & IPmax ) I 



(4.43) 



and correspondingly for p min , where a is the function from eq. (4.35). 

To summarize, all spectral functions corresponding to real corrections have 2-dimensional 
integral representations: 



Prl 



tt/CV^o) 



(4tt) 4 £; 
sign(a) 

F(po,q) 



dp Q dq n (k - p ) n (q) n (p - q) 



In 



Po 



p 2 . 

rmin 



Po 



Pi 



+ In 



[F(po,q) + |p max a|] 2 - (q - k ) 2 (p 2 - ptjf 



[F(po, 9) + |p min a|] 2 - (g - feo) 2 (Po - Pmin) 2 



(4.44) 



and correspondingly for p r2 , p r3 and p r4 . It can be recalled from eq. (4.28) that the other 
channels come with an overall minus sign, and we also find that sign(a) = — 1 for the channels 
(r2) and (r4). The integration range fi rl + Q r2 + ^r3 + ^r4 i s as given in fig. 4. 

Although well suited for numerical handling, the integral representation in eq. (4.44) re- 
mains fairly complicated in practice as long as A 7^ 0. For A — > 0, the expression simplifies 
but is also logarithmically divergent. However, when we sum the result together with vir- 
tual corrections, to which we now turn, the divergences cancel; the resulting expressions are 
presented in sec. 6. 



4 The function a can be positive or negative but not change its sign within the ranges considered. 
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5. Virtual corrections 

For ko > k > 0, the virtual corrections are contained within the terms denoted by (vl)— 
(v3) in eq. (3.2), as well as in their reflections (02 -H- cri,<74 o 03). Each of the terms is 
factorized into two structures. In the first one the angular integration is immediately doable, 
and thereby we can simplify the radial integration measure into 



27r5(-k + e_ + e 



pk) 



ie p € pk 



1 + n a ,(e pk ) + n (e ) 



n ao(ko) 
8-rrk 



dpn (k -p)n (p) . (5.1) 



As far as the q-integrals are concerned, their angular parts can be carried out with the help 
of a Feynman parameter, s £ (0, l). 5 Considering first the case (vl), we note that 



s[(Q - Vf - A 2 ] + (1 - s)(Q - K,) 2 = m 2 + 2q • e s 



(5.2) 



where m 2 = (1 — s)K? — s\ 2 — 2qo[sp + (1 — s)ko] and e s = sp + (1 — s)k. Denoting by z the 
angle between q and e s , we trivially get 



+1 



dz 



1 



1 



_i 2 (m 2 + 2qe s z) 2 m\ — Aq 2 e 2 s 



(5.3) 



Furthermore, by making use of the constraint ka = p + e k implied by eq. (5.1), the length of 
e s is given by e 2 = [sp + (1 — s)ko] 2 — (1 — s)/C 2 . The denominator in eq. (5.3) is thus a 2nd 
order polynomial in s, and the integral over s is also doable: 



ds 



1 



In 



P- V^ 



up 



^s 2 -2us + p 2a/ 'v 2 - pp v-p+^/v 



pp 



As a final step, the terms with go = ± e q can formally be combined by making use of 



1 

— + TL r 

2 ° 



L 2 + n ^ ( i) 



(5.4) 



(5.5) 



Thereby the term (vl) can be cast in the form 

dpn a (k -p)n a (p) 



Pvl 



*JC 4 n- l (ko) r*+ 



(4vr) 4 /c 



A;_ 



dq 



+ n a Aq) 



1 



7"(p,i 



In 



2/C 2 g(g - p) - A 2 (/C 2 - 2fco<7) + 2 9 J-(p, ?) 



2/C 2 g(g - p) - \ 2 {K, 2 - 2k q) - 2qT(p, 



(5.6) 



The vacuum part, which is represented by the factors | in eq. (3.2), could be integrated explicitly by 
making use of Lorentz invariance (the result is shown in eq. (5.10) below), however for us it is convenient 
to handle it together with the thermal contributions. For a number of other master spectral functions the 
vacuum part is divergent at large |q|, but it turns out that even in those cases it is convenient to handle it 
together with the thermal contributions for moderate |q| <fco; only the asymptotics at |q| ^> ko needs to be 
handled separately with a proper ultraviolet regularization. 
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where T is the same function as appeared in the real corrections, given by eq. (4.42). The 
argument of the square root in T is always positive, so that the integrand is well-defined; 6 
however, the vacuum part of the g-integral in eq. (5.6) is divergent on its own, and the integral 
should only be carried out for the sum of the three terms (vl)-(v3). 

The second structure (v2) can be reduced to the first one by substituting variables as 
p — )■ k — p, q — t- k — q in eq. (3.2). Effectively, this corresponds to an interchange 04 -H- u\ 
in eq. (5.6). However, we can subsequently also substitute p — >■ ko — p, q — > ko — q, and given 
that J-{p, q) of eq. (4.42) is invariant in this transformation, the result reads 



Pv2 



7r/C 4 n^(fc ) 



1 



(4vr) 4 £: j k _ 

2JC 2 (q-k )(q-p) 



In 



d P n a 4 ( k o-p)n ai (p) / dg 
J — 00 

\ 2 [JC 2 + 2k (q 



[ \ + "<* 



ko)] + 2(k -q)J r (p,q) 



F(p, q) 

For the third term 
with 

s(V-Q 

where now m 2 = X 2 — 2qo(p 
go as before, noting that e 2 



(5.7) 



2K?(q - ko)(q -p)~ A 2 [/C 2 + 2k (q - ko)] - 2(k - q)T(p,q) 

v3), we can substitute q — > p — q in eq. (3.2). Eq. (5.2) gets replaced 



/C) 2 + (l-s)(V-Q) 2 =m 2 + 2q-e s 



(5.8) 



sko 



p — sk, and qo = ±E„. The subsequent steps 



p + s(JC — 2pko) + s k . Afterwards, it is convenient to 



again return to the "original" variables; this can be implemented by first taking E q as an 
integration variable instead of q, and then introducing a "new" q as E q = \p — q\. In this way 
the result can be cast in a form reminiscent of eqs. (5.6), (5.7): 



Pv3 



7r/C 4 n- ( ,(fco) 
(4vr) 4 £, 

IC 2 (q 



dpn (k -p)n a (p) 



k- 



p-X 



+ 



F(p, l) 



In 



P+X. 
X 4 



dq 



\ +n ^ 



(p-q) 



p)2 _ x 2 [K 2 + 2qp -k (q+p)}-^- ^f(q-p) 2 -\ 2 F(p, 



K 2 (q - p) 2 - A 2 [/C 2 + 2qp - k (q + p)] - f + y/(q - p) 2 - \ 2 F(p, 



(5.9) 



The expressions in eqs. (5.6), (5.7), (5.9) contain no expansion with respect to A. We have 
crosschecked numerically in the small-A regime that picking up the vacuum factors (the "g "s 
if the n's are expressed with positive arguments), their sum is integrable and reproduces 
within the achieved resolution a known Lorentz-invariant integral: 

1 



IA 



Q \[(Q-P) 2 + \ 2 ]Q 2 (Q 



1 



(4vr) 2 /C 2 



7T 

~6 



ln z 



K) 2 
\ 2 



K 2 + A 2 



P 2 =0,K 2 

- LL 



-K 2 ,K-P~- 
X 2 

/C 2 + A 2 



2 



(5.10) 



6 It may be noted, however, that the integrand has a non-trivial structure at small q, with the argument 
of the logarithm having a zero at q = — A 2 /(4p). This (integrable) singularity cancels against a corresponding 
one from real corrections, related to structure near the boundaries between the regimes (rl) and (r3) in fig. 4, 
cf. footnote 7. 
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6. Cancellation of divergences 

The results in sees. 4 and 5 are both divergent if we attempt to send A — > (cf. eq. (5.10) for 
the vacuum part). Their sum, however, remains finite as we now demonstrate. 

The expressions obtained have two kinds of "singularities" . In terms of fig. 4, with variables 
suitably renamed for virtual corrections, the integrands have a non-trivial structure around 
the boundaries separating the different channels, i.e. q = po, q = 0, and po = k$. In addition, 
away from the boundaries, the integrands in general diverge as A — > 0. 

The strategy we adopt is to stay away from the boundaries, for instance by setting a band 
of width S around them, and taking the limit A — > within the domains. 7 We then verify 
the cancellation of the corresponding divergences within the domains, and that the resulting 
integrand remains integrable even after ultimately setting <5 — >■ 0. 

Turning first to the real corrections, eq. (4.44), we simplify the notation from now on by 
renaming 

Po^p. (6.1) 

Then the following expressions are obtained for the integrand of eq. (4.44) inside the domains 
of fig. 4: 

1 



sign(a) 

F{p, q) 



A^O 



(a) 



-(1): 
(b): 

(b): 

(h) = (h) = (j) : 

(c): 

(d): 

(e) = -(f): 



(p-q)JC< 

1 
(p-q)& 

1 
(p-q)& 

1 
(p-q)& 

1 
(p-q)IC 2 

1 

(p-q)IC 2 



In 



K\p-qf 



(p-q)& 



In 



In 



In 



In 



In 



In 



X 4 pq(p - k )(q - k 
K,\p-qf(q-k + ) 



\ 2 q{q 

ic 2 (p- 



k )(p 

q) 2 {p- 



-k, 
k- 



X 2 p(p - k )(q - 
(P ~ k-){q - k+ 



fc_ 



(p 



k+)(q-k.) 

K\p 



qf 



\ 4 p(q-k )(p- 

q(p - k ) 



k-)(q-k+) 



(P 



-k+)(q-k-) 

lO{p-qf(p 



\ 2 (q - k )(p - k-)(p - k + ) 



(6.2) 
(6.3) 
(6.4) 
(6.5) 
(6.6) 
(6.7) 
(6.8) 



7 Actually this requires a somewhat more careful justification, given that according to eq. (4.28) there is 
a phase space distribution at each boundary which would diverge if it were bosonic. It turns out that if 
the integrand of eq. (4.44) is evaluated within the domains adjacent to the boundaries and we make the 
substitution in eq. (6.1), then it cancels exactly against an integrand in one of the virtual corrections, namely 
that containing the same potentially divergent phase space distribution. In other words, the 2nd row of 
eq. (4.44) evaluated within the domains (a) and (1) of fig. 4 cancels against the 2nd row of eq. (5.9); within 
(e) and (f) against that of eq. (5.7) once the latter is reflected in p <-> q; and within (e) and (f) against that 
of eq. (5.6). Approximate forms of these cancellations can be seen by contrasting eqs. (6.2) and (6.15); (6.8) 
and (6.14) [with p <-> q]; as well as (6.9) and (6.13). 
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Figure 5: Integration ranges in the (p, q)-p\ane for virtual corrections. The delimiting curves are 
p — k—, p — k+, q = k-, and q = k+. Narrow bands have been cut off around the lines q = p, q = 
and p = ko- 



(e) = -(f) 
(g) = -(g) 
-(i) = ~(k) 
■Q) = "(k) 





1 


(p- 


-g)/C2 
1 


(p- 


-g)/C2 
1 


(p- 


-g)/C2 
1 



(p-g)/C 2 



/C 2 (p — q) 2 q 



X 2 p(q - k-)(q - k+) 
q(p - k ) 



(p-k-)(q-k+) 

JC 2 (p-q) 2 (q-k_) 



X 2 q(q - k )(p - k-) 
K 2 (p-q) 2 {p-k + ) 



X 2 p(p - k Q )(q - k+) 



(6.9) 
(6.10) 
(6.11) 
(6.12) 



These are multiplied by phase space distributions as indicated by eq. (4.28), and in addition 
the symmetrization (02 f* <J\ , 0"4 O 173) from the last line of eq. (3.2) needs to be included. 
Similar expressions are obtained for the virtual corrections, eqs. (5.6), (5.7), (5.9): 



1 



F(p, q) 



in 



A^O 



(vi: 



(v2^ 



(v3) 



1 



X 2 p(q - fc_)(g - k+) 



(p- 


- q)IC 2 " 
- 1 In 


(p- 


- q)JC 2 
1 In 


\p- 


- q\IC 2 



K?q(p — q) 2 
X 2 (p-k )(q-k_)(q-k + ) 



K?(q-ko)(p-q)< 
X A pq(p- ko)(q- k ) 



IC 4 (P 



(6.13) 
(6.14) 
(6.15) 



Again the symmetrization (02 f* 01, 04 -O- 03) needs to be included. 

In order to combine eqs. (6.13)-(6.15) with the real corrections, it is beneficial to rename 
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integration variables. By making use of 

n a (-e) = -1 - n a (e) , n a (e)n T (5 - e) = n ar {5) [l + n a (e) + n T {5 - e)] , (6.16) 

and recalling the identities o"2<75 = a±, cr^a^, = 04, following from eq. (2.6), we can write 

2 + n a 2 (Q) = n ai(p) n a 2 (Q) n a 5 (P ~ l) ~ ^ ~ n »B (P ~ ?) ' ( 6 - 17 ) 

- + n CT3 (/c - ?) = n~l{ko ~ P)n a3 (k - q)n as (q - p) + - + n CTS (p - g) . (6.18) 

Furthermore, in the first term on the right-hand side of eq. (6.18), in which n ff3 (ko — q) 
appears, we exchange variables as p -H- q. In the terms from (02 f* o\ , 04 ■<->■ 03), we can also 
do this in the terms involving ^ + n CTr (p — q) for more symmetry. 

The effect of these rewritings is that logarithms of A disappear from terms involving ^ + 
n a 5 (P ~ l) [essentially its coefficient is given by —(6.13) + (6.14) + (6.15)]. Logarithms of 
A do not cancel from the coefficients of the first terms of eqs. (6.17), (6.18) but, combining 
with the other phase space distributions from eqs. (5.6), (5.7), are seen to come with the 
same "weight functions" as in the real corrections. To be explicit, the virtual corrections 
[p v = p vl + p v2 + p v3 + (o"2 ** o\, 04 ** 03)] can be represented within the domains shown in 
fig. 5 as 



(47T) 4 fc / > v (/C) 



dpdg 



n U4 (ko - p)n a2 (q) + n (k - p)n ffl (q) n a . (p - q) 



X 2 p(q - k-)(q - k+) 



(p~q)IC 2 



+ n * 4 (k - p)n (p) - + n (p - q) 




+ 



dpdq 



(p- q )W (q-k.)(q-k + ) 



n ai (k - p)n a2 (q) + n (k - p)n ffl (q) n (p - q) 



\ 2 (q-ko)(p-k^){p-k + ) 




+ n a, ( k ~ l) n ao (?) o + U o. (P ~ ?) 



(p-q)K, 2 



(6.19) 



Combining eq. (6.19) with eqs. (6.2)— (6.12), the latter multiplied by phase space distribu- 
tions according to eq. (4.28), and adding for the real corrections terms from the symmetriza- 
tion (o"2 <-> <7i , 04 (73), all logarithms of A are now seen to cancel. The resulting integrand 
is sufficiently well-behaved around the boundaries to be integrable [cf. eq. (7.9) and sec. 8]. 
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Figure 6: Integration ranges in the (p, g)-plane for the final answer, eq. (7.9). The delimiting curves 

are p = k-, p = k+, q = p, q = p + k-, q = p + k+ , q = ko — p, as well as q = fc+ . 

7. Final result 

In order to collect together a final result, we remove redundant symmetries from the expres- 
sion. It is suggested already by eq. (4.42) that the two substitutions 



p-H- k -q; p^ q 



(7.1) 



corresponding to reflections across the axes q = k^—p and q = p, respectively, may be helpful 
in this respect. Indeed, the first of these symmetries is manifest in the real corrections, 
and can consequently be seen to transform the "coefficient functions", eqs. (6.2)-(6.12), into 
each other. The phase space distributions are in general not invariant but now become 
symmetrized with respect to their indices (concretely a symmetry in (<7i o <74,o"2 ^ 03) 
which was so far hidden becomes explicit). 

A symmetry in p f> g is not "inherent" to the expressions. It is useful to make this 
reflection, however, because it explicitly regulates principal value integrations across p—q = 0. 
After these reflections, the integration range is as shown in fig. 6. 

To present a final result, we undertake one more rewriting of the phase space distributions. 
By making use of eqs. (5.5), (6.16), the structures appearing in the real corrections are 
expressed as 

n ai {k - p)n a2 (q)n a5 (p - q) = n a4 (k - p)n ai {p) { -+n a2 (q) - -+n U5 (q-p) }, (7.2) 
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where n ar has been taken with a positive argument. We then introduce the "weight functions" 



Wl = n a4 (k - p)n ai (p) - + n a5 (q - p) + (a u a 4 ,a 5 ) + (a 2 ,a :i ,a 5 ) + (a 3 ,a 2 ,a 5 ) 

^2 = n ff4 (k - q)n ai (q) - + n ff5 (q - p) + (<n ,<t 4 ,<t5) + (o- 2 ,<r 3 ,o- B ) + (0-3,0-2 ,* B ) , 
w 3 = n (k -p)n (p) - + n (q) + (0-1,0-4,03) + (o 2 ,o 3 ,o 4 ) + (o- 3 ,o- 2 ,oi) , 



iO.i 



.2 

n a4 (k - q)n ai (q) - + n a2 (p) + (0-1,0-4,0-3) + (o- 2 ,o 3 ,o 4 ) + (o 3 ,o 2 ,oi) 

To characterize the "coefficient functions" it is convenient to define the ratios 

q — k_ _ Q — k+ _ q 



X- 



7T_ = 



q-k 

p — k- 
p-k 



X+ 



q-k 

P-K 
p-k 



Xo 



> ^0 



q-k Q 

P 
p-k 



(7.3) 
(7.4) 
(7.5) 
(7.6) 

(7.7) 
(7.8) 



Then the full result becomes 

(47T) 4 fcp x (£) 



7r/C 2 n CT0 (ko) 



+ 



+ 



+ 



+ 



+ 



+ 



+ 



+ 



+ 



dpdqf 

< (wi - w 2 )m 

n, , Q ~ V I 



7T IT, 



+ W3I1I 



Xo 



k+b 



dpdgf 

< wi In 



P 



X^+ 



X-X\ 



dpdqf, . 

<^ (wi + W2 - W3 + W4j In 



X-X+ 
+ (W2 + W4) In 

X-7T+ 



X-X+7T-7T+ 

X-7T+ 



7T n 



+ W3ln 



+ a^ln 

Xo 



X-X+TT-^ 



7Tn 



X-7T+ 






dpd^f 

■<wiln 



c 

dpdq ( 



7T_7T 1 



X-X+ 



Wi — W3) In 



+ W2I11 

X-7T+ 



X+7T- 

^oX 



7r_7Tj 



+ w 3 ln 



Xo 



X-X+7T-7T+ 



+ u^ln 



X+7T- 



7Tn 



k+c+f 



dpdq( 

< uj\ In 



Xo 



7T. 



X-X+ 



X+7T- 

+ W3 In 



+ (u> 2 + W4) In 

Xo 



7T, 







X+7T- 



X-7T+ 



dpd^f 

^ (cj 2 + W4) In 



n g q~P 

dpdq ( 
<X Q-P 



vr, 







I 



(oj 2 + W4) In 



X-7T+ 
X+7T- 



dpdq( 

< wi in 



P 



k+i 
dpdq ( 



Xo 



7T, 



X-X+ 



X-7T+ 

+ (W2 + 0J4) In 



X+7T- 



7T, 







+ u; 3 ln 



Xo 



X-7T+ 



n.Q-P 



\ 



[ui -w 3 )ln 



X+7I"- 



(7.9) 
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Figure 7: Left: A comparison of the spectral function at non-zero momentum (k/T — 0.1, 0.3, 1, 
3, 4, 6, 8 as indicated in the figure) with the zero-momentum limit from eq. (A. 57) of ref. [20], for 

(crio-4£r 5 ) = (+ + +). For the case (cria^CTs) = (H h) only non-zero momenta are shown. In both 

cases we have restricted to ko > k + 0.001T. Right: A magnification of the case (ctio^o^) = (+ + +). 

where the Q's denote different domains as labelled in fig. 6. Note that only one among the 
ranges Slj+ C and Oi + g gets realized at a time. 



8. Numerical evaluation 

The expression in eq. (7.9) is finite and can be evaluated numerically. In fig. 7 the outcome 
is compared with its limiting value at k = 0, determined for the case (<7i 04(75) = (+ + +) 
in ref. [20]. The results are seen to agree for k$ 3> k; in fact, even when this inequality is 
not satisfied, the zero-momentum limit yields a surprisingly good approximation of the full 
result. 

It is important to realize that the spectral function p x is in general non-trivial in the 

j 

vicinity of the light cone, cf. fig. 7, rather than vanishing as ~ K, 2 as one could expect from 
dimensional reasons at zero temperature. This fact leads ultimately to the breakdown of the 
loop expansion for /C 2 <(gT) 2 [2-loop diagrams may give a larger contribution than 1-loop 
ones despite the overall suppression by g 2 ], and leads to the need to resum a set of diagrams 
for obtaining the correct result in the ultrarelativistic regime. 

Another comparison can be made with the non-relativistic asymptotics, determined in 
ref. [16]. The leading term, proportional to /C 2 , corresponds to the zero-temperature limit 
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Figure 8: Left: Spectral functions for various statistics (cio^o^), as a function of M 2 = K, 2 , with 
k = k &v determined from eq. (8.3). The results are compared with two orders of OPE asymptotics 
from eq. (8.1). Right: The ratio [exact - L0 OPE ]/NL0 OPE , showing that numerical results agree with 
the asymptotic ones from eq. (8.1) for M >30T. 



and vanishes for the spectral function in question. The two first non-zero terms read 

>4 _ , „„ r^ 4 



Pi- 



■—( 

16vr J p 



Eti^ + 2n. 



<Jr, 
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i=X 
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P P ~ 12 ' J p P " 24 ' J/"' 
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In fig. 8 the two orders shown in eq. (8.1) are referred to as "L0 OPE " and "NL0 OPE ", respec- 
tively. In order to carry out the comparison in a somewhat realistic setting, we introduce a 
phenomenological average momentum through 



Vk 2 +M 2 



kl(M) 



f™dkk 4 eM- T ) _ 3MTK 3 (f) 
f^dkk^ex V (-^±Ml)~ K 2 (f) 



(8.3) 



This should be understood just as a rough guideline; for instance we have employed a Boltz- 
mann weight so that the same value can be used for any statistics. In any case, as can clearly 
be seen in fig. 8(right), the correct limits are reached for all statistics considered, if only quite 
deep in the non-relativistic regime. 
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9. Summary and outlook 

The purpose of this paper has been to suggest a general strategy for determining 2-loop 
thermal spectral functions at non-zero energy and momentum in the rest frame of a heat bath. 
As has been demonstrated with the example of the most complicated "master" structure, the 
result can be reduced to a convergent 2-dimensional integral, eq. (7.9), within a domain 
shown in fig. 6. For other master spectral functions, the domain remains the same but the 
"coefficient functions" change; in addition, due to less symmetry, the "weight functions" split 
up into a larger set of independent ones (the weight functions are defined as polynomials of 
the phase space distributions incorporating all the temperature dependence). 8 

The specific master spectral function studied, defined by eqs. (2.5), (3.1), is peculiar in that 
it vanishes in the zero-temperature limit (this comes about through a complete cancellation 
of real and virtual corrections, and is reproduced by our numerical results). Therefore it 
is natural to express it as p x = T 2 cp(k + /T,k^/T), where k± = (fco ± k)/2 and is a 
dimensionless function. The function <f> has a finite limiting value in the non-relativistic limit 
fc + ,fc_ 3> ttT, but a non-trivial structure in the relativistic regime &+,&_ ~ irT, cf. fig. 7. 
(In the regime k- <C irT the naive loop expansion of thermal field theory breaks down and 
needs to be resummed through effective field theory techniques.) 

Several extensions of the current investigation can be envisaged. The most obvious chal- 
lenges are to work out similar results for the other master structures appearing in eqs. (2.4), 
(2.9) and then to compile results for the physical observables discussed in sec. 2. It might 
also be interesting to extend the results to a situation where some of the propagators are 
massive; this would be relevant for the cosmological applications reviewed in ref. [14]. (As 
has been demonstrated with a particular non-zero mass here, it may be possible to reduce 
the result to a 2-dimensional integral even in the presence of masses.) Perhaps it would be 
nice to understand analytically the behaviour in the regime K? <C (ttT) 2 . In addition the 
question could be posed whether, possibly with the price of introducing one further integra- 
tion variable, the final result in eq. (7.9) could be cast in a more compact and transparent 
form. Last but not least, the computation presented involved a fair amount of error-prone 
hand work, so that an independent crosscheck, perhaps involving other integration variables 
and/or another intermediate infrared regulator, would be more than welcome. 
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